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Abstract. We study numerically solutions to the Korteweg-de Vries and 
Camassa-Holm equation close to the breakup of the corresponding solution 
to the dispersionless equation. The solutions are compared with the properly 
rescaled numerical solution to a fourth order ordinary differential equation, 
the second member of the Painleve I hierarchy. It is shown that this solu- 
tion gives a valid asymptotic description of the solutions close to breakup. We 
present a detailed analysis of the situation and compare the Korteweg-de Vries 
solution quantitatively with asymptotic solutions obtained via the solution of 
the Hopf and the Whitham equations. We give a qualitative analysis for the 
Camassa-Holm equation 



1. Introduction 

It is well known that the solution of the Cauchy problem for the Hopf equation 

(1.1) u t + 6uu x = 0, u(x,t = 0) =u (x), xeR, t e R+ 

reaches a point of gradient catastrophe in a finite time. The solution of the vis- 
cosity or conservative regularization of the above hyperbolic equation display a 
considerably different behavior. Equation (|1.1[) admits an Hamiltonian structure 

u t + {u(x),H a } = u t + d x — r^r = 0, 

0U(X) 

with Hamiltonian and Poisson bracket 

H Q = u 3 dx, {u(x),u(y)} = 5'(x-y), 
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respectively. All the Hamiltonian perturbations up to the order e of the hyperbolic 
equation (11. ip have been classified in |10j . They are parametrized by two arbitrary 
functions c(u), p(u) 

e 2 

lit + 6u u x + — [2cu xxx + 4c'u x u xx + c"u x ] + e 4 [2pu xxxxx 

(1.2) 

+2p' (bu xx u xxx + 3u x u xxxx ) + p"{7u x u 2 xx + 6u 2 x u xxx ) + 2p"'u 3 x u xx \ = 0, 

where the prime denotes the derivative with respect to u. The corresponding Hamil- 
tonian takes the form 



H 



24 



<? x + e A p{u)u 2 x:i 



dx 



For c(u) = 12, p(u) = one obtains the Korteweg - de Vries (KdV) equation 
u t + 6uu x + e u xxx = 0, and for c(it) = 48it and p(u) = 2u the Camassa-Holm 
equation up to order e 4 ; for generic choices of the functions c(u), p(u) equation 
(jl.2p is apparently not an integrable PDE. However it admits an infinite family of 
commuting Hamiltonians up to order 0(e 6 ). 

The case of small viscosity perturbations of one-component hyperbolic equa- 
tions has been well studied and understood (see pQ and references therein), while 
the behavior of solutions to the conservative perturbation (|1.2p to the best of our 
knowledge has not been investigated after the point of gradient catastrophe of the 
unperturbed equation except for the KdV case, |18l 1231 [7]- 

In a previous paper [13] (henceforth referred to as I) we have presented a 
quantitative numerical comparison of the solution of the Cauchy problem for KdV 

(1.3) u t + 6uu x + e 2 u xxx = 0, u(x, 0) = u (x), 

in the small dispersion limit e — > 0, and the asymptotic formula obtained in the 
works of Lax and Levermore |18j . Venakides [23] and Deift, Venakides and Zhou [7] 
which describes the solution of the above Cauchy problem at the leading order as 
e — > 0. The asymptotic description of |18j .[7 ] gives in general a good approximation 
of the KdV solution, but is less satisfactory near the point of gradient catastrophe 
of the hyperbolic equation. This problem has been addressed by Dubrovin in [10] , 
where, following the universality results obtained in the context of matrix models 
by Deift et all [8], he formulated the universality conjecture about the behavior of 
a generic solution to the Hamiltonian perturbation (|1.2p of the hyperbolic equation 
(|l.ip near the point (x c ,t c ,u c ) of gradient catastrophe for the solution of (jl.ip . He 
argued that, up to shifts, Galilean transformations and rescalings, this behavior es- 
sentially depends neither on the choice of solution nor on the choice of the equation. 
Moreover, the solution near the point (x c ,t c ,u c ) is given by 
(1.4) 

u(x, t,e)~u c + a e 2/7 U (b e- 6/7 (x -x c - 6u c (t - i c )); ccT^\t - t c j) + O (e 4/7 ) 

where a, b, c are some constants that depend on the choice of the equation and the 
solution and U = U(X;T) is the unique real smooth solution to the fourth order 
ODE 



(1.5) X = 6TU 



U 3 + {\u 2 x + UU X x) + ^U xxxx 
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which is the second member of the Painleve I hierarchy. We will call this equation 
PI2. The relevant solution is characterized by the asymptotic behavior 

(1.6) U(X,T) = T{X)i T^r + OiX- 1 ), X -> ±00, 

X 3 

for each fixed T £ R. The existence of a smooth solution of (|1.5p for all X, T G R 
satisfying (|1.6p has been recently proved by Claeys and Vanlessen [4] . Furthermore 
they study in [5] the double scaling limit for the matrix model with the multicritical 
index and showed that the limiting eigenvalues correlation kernel is obtained from 
the particular solution of (11. 5p satisfying (|1.6p . This result was conjectured in the 
work of Brezin, Marinari and Parisi [2] . 

In this paper we address numerically the validity of (| 1 .4[) for the KdV equation, 
and we identify the region where this solution provides a better description than the 
Lax-Levermore, and Deift-Venakides-Zhou theory. As an outlook for the validity of 
(|1.4p for other equations in the family (|1.2p , we present a numerical analysis of the 
Camassa-Holm equation near the breakup point. While the validity of (|1.4p can 
be theoretically proved using a Riemann-Hilbert approach to the small dispersion 
limit of the KdV equation [7] and recent results in [§], [4], [5], for the Camassa-Holm 
equation and also for the general Hamiltonian perturbation to the hyperbolic equa- 
tion (jl.ip . the problem is completely open. Furthermore for the general equation 
(|1.2p . the existence of a smooth solution for a short time has not been established 
yet. An equivalent analysis should also be performed for Hamiltonian perturbation 
of elliptic systems, in particular for the semiclassical limit of the focusing nonlinear 
Schrodinger equation |16j.|21j. 

The paper is organized as follows. In section 2 we give a brief summary of 
the Lax-Levermore, and Deift-Venakides-Zhou theory and the multiscale expansion 
(|1.4p . In section 3 we present the numerical comparison between the asymptotic 
description based on the Hopf and Whitham solutions and the multiscale solutions 
with the KdV solution. In section 4 we study the same situation for the Camassa- 
Holm equation. In the appendix we briefly outline the used numerical approaches. 



2. Asymptotic and multiscale solutions 

Following the work of [18] , [23j and [7] , the rigorous theoretical description of 
the small dispersion limit of the KdV equation is the following: Let u(x,t) be the 
zero dispersion limit of u(x,t, e), namely 

(2.1) u(x, t) = lim u(x, t, e). 

e— >0 

1) for < t < t c , where t c is a critical time, the solution u(x, t, e) of the KdV 
Cauchy problem is approximated, for small e, by the limit u(x,t) which solves the 
Hopf equation 

(2.2) u t + 6uu x = 0. 

Here t c is the time when the first point of gradient catastrophe appears in the 
solution 



(2.3) 



u(x,t)=u (Q, x = 6tu (0 + £, 
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of the Hopf equation. From the above, the time t c of gradient catastrophe can be 
evaluated from the relation 

1 



tr 



2) After the time of gradient catastrophe, the solution of the KdV equation is 
characterized by the appearance of an interval of rapid modulated oscillations. 
According to the Lax-Levermore theory, the interval [x~(t), x + (t)] of the oscillatory 
zone is independent of e. Here x~(t) and x + (t) are determined from the initial 
data and satisfy the condition x~(t c ) = x + (t c ) = x c where x c is the x-coordinatc 
of the point of gradient catastrophe of the Hopf solution. Outside the interval 
[x~~ (t), x + (t)] the leading order asymptotics of u(x, t, e) as e — * is described by the 
solution of the Hopf equation (|2.3p . Inside the interval [x~ (t) , x + (t)] the solution 
u(x,t, e) is approximately described, for small e, by the elliptic solution of KdV 



(2.4) u(x, t, e) ~ G + 2e 2 ^ log 9 {^E^[ x - 2i(ft + ft + ft) - q] ; T 

where now u = u{x, t) takes the form 

(2.5) u = ft + ft + ft + 2a, 

fort a r t - K '^ 2 

(2.6) a = -ft + (ft -W^H, T = i-f-rr-r, s z = — 

K(s) K(s) ft - ft 

with K(s) and E(s) the complete elliptic integrals of the first and second kind, 
K'(s) — K(yl — s 2 ); 6 is the Jacobi elliptic theta function defined by the Fourier 
series 

d{z;T) = ^ e ™ 2 ^+2««^ 

riGZ 

For constant values of the ft the formula (|2.4p is an exact solution of KdV well 
known in the theory of finite gap integration [15] , [9j . However in the description 
of the leading order asymptotics of u(x, t, e) as e — > 0, the quantities ft depend on 
x and t and evolve according to the Whitham equations [24) 

^ft + ^^ft = 0, i= 1,2,3, 
where the speeds Vi are given by the formula 

(2.7) Vt= J l ^f;-^ + 2(ft+ft + ft), 

ft + a 



with a as in (j2.6j) . Lax and Levermore first derived, in the oscillatory zone, the 
expression (|2.5p for u = u(x, t) which clearly does not satisfy the Hopf equation. 
The theta function formula (j2.4[) for the leading order asymptotics of u(x,t,e) as 
e — * 0, was obtained in the work of Venakides and the phase q — ?(ft, ft, ft) was 
derived in the work of Deift, Venakides and Zhou [7J, using the steepest descent 
method for oscillatory Riemann-Hilbert problems 6 

(2.8) g(ft,ft,ft) = — 7= / / fljUW- " " 



where f—(y) is the inverse function of the decreasing part of the initial data. The 
above formula holds till some time T > t c (see [7J or I for times t > T). 
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3) Fei-Ran Tian proved that the description in 1) and 2) is generic for some time 
after the time t c of gradient catastrophe |20j . 

In I we discussed the case uq(x) = — sech x in detail as an example. The main 
results were that the asymptotic description is of the order 0(e) close to the center of 
the Whitham zone, but that the approach gives considerably less satisfactory results 
near the edges of the Whitham zone and close to the breakup of the corresponding 
solution to the Hopf equation. In the present paper we address the behavior near 
the point of gradient catastrophe of the Hopf solution in more detail. In Fig. Q] we 
show the KdV solution and the corresponding asymptotic solution as given above 
for several values of the time near the critical time t c . It can be seen that there are 
oscillations before t c , and that the solution in the Whitham zone provides only a 
crude approximation of the KdV solution for small t — t c . The situation does not 
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Figure 1 . The blue line is the solution of the KdV equation for the 
initial data uq(x) — —1/ cosh 2 x and e = 10~ 2 , and the purple line 
is the corresponding leading order asymptotics given by formulas 
(|2.3p and (|2.4p . The plots are given for different times near the 
point of gradient catastrophe (x c ,t c ) of the Hopf solution. Here 
x r ~ -1.524, t r ~ 0.216. 



change in principle if we consider smaller values of e as can be seen from Fig. [5] 
The solution shows the same qualitative behavior as in Fig. [Tl just on smaller scales 
in t and x. 

2.1. Multiscale expansion. We give a brief summary of the results in |10j 
relevant for the KdV case we are interested in here. Near the point of gradient 
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FIGURE 2. KdV solution and asymptotic solution for e = 10 3 
close to the breakup time. 

catastrophe (x c , t c ,u c ), the Hopf solution is generically given in lowest order by the 
cubic 

(2.9) x-x c ~ 6(t-t c )u- k(u-u c ) 3 , k = -f"(u c )/6, 

because 6t c + f'_(u c ) — and f"(u c ) = 0. Here f-(u) is the inverse of the decreasing 

SHh 

part of the initial data uq(x). Now let us consider hk = — j — where Hk are the 

ou 

KdV Hamiltonians such that h k = u k+2 /{k + 2)! + 0(e 2 ). We have 

u 2 e 2 1 e 2 e 4 

h-i=u, h = — + —u XXl hi = -(u 3 + ~^{u x + 2uu xx ) + — u xxxx ), 

and the KdV equation is obtained from u t + 6d x ho = 0. Then 

x = 6tu + aoho + a±hi + . . . a,khk, 

is a symmetry of the KdV equation [11] . Setting ao = 0, a\ = —f"'(u c )/6 = k 
and Ofc>2 = 0, and making the shift t — > t — t c , u — > u — u c and the Galilean 
transformation x —> x — x c — 6(i — t c )u c we arrive at the fourth order equation of 
Painleve type 
(2.10) 



x— x c — 6{t— t c )u c — 6(t~t c )(u—u c ) — k 



u 2 



(u - u c ) + e + (u — u c )u xx ) + —u xxxx 



which is an exact solution of the KdV equation and can be considered as a pertur- 
bation of the Hopf solution (|2.9I) near the point of gradient catastrophe (x c , t c , u c ). 
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The solution u(x,t,e) of ([2"TU|> is related to the solution U(X,T) of (|L5|) by the 
rescalings 

(2.11) U (x,i,e) = Mc +(^) 2/ V(X,T) 
where 

£ — X r — 6u r (t —t r ) „ i — fr 

(2.12) X= t—r-P cl , T = ^-. 

According to the conjecture in [10] . the solution (|2.11[) is an approximation modulo 
terms 0(e 7 ) to the solution of the Cauchy problem (|1.3|> for (x, t, u) near the point 
of gradient catastrophe (x c ,t c ,u c ) of the hyperbolic equation (12.21) . 

3. Numerical comparison 

In this section we will present a comparison of numerical solutions to the 
KdV equation and asymptotic solutions arising from solutions to the Hopf and 
the Whitham equations as well as the Painleve 12 equation as given above. Since 
we control the accuracy of the used numerical solutions, see I, [17] and the ap- 
pendix, we ensure that the presented differences are entirely due to the analytical 
description and not due to numerical artifacts. We study the e-dependence of these 
differences by linear regression analysis. This will be done for nine values of e be- 
tween 10 _1 and 10~ 3 . Obviously the numerical results are only valid for this range 
of parameters, but it is interesting to note the high statistical correlation of the 
scalings we observe. We consider the initial data 

uo(x) = — 1/ cosh 2 a;. 

For this initial data 

(3.1) x c = -^+log((V3-l)A/2), t c = ^, u c =-2/3. 

3.1. Hopf solution. We will first check whether the rescalings of the coor- 
dinates given in ([2.11|1 arc consistent with the numerical results. It is known that 
the Hopf solution provides for times t <C t c an asymptotic description of the KdV 
solution up to an error of the order e 2 . This means that the Loo-norm of the differ- 
ence between the two solutions decreases as e 2 for e — > 0. For t — 0.1 we actually 
observe this dependence. More precisely this difference can be fitted with a 
straight line by a standard linear regression analysis, — log 10 A M = — alog 10 e + b 
with a = 1.9979, with a correlation coefficient of r = 0.99999 and standard error 
a a = 4.1 * IQ- 3 . 

Near the critical time t c this picture is known to change considerably. Dubrovin's 
conjecture |10| presented above suggests that the difference between Hopf and KdV 
solution near the critical point should scale roughly as e 2 / 7 . In the following we 
will always compare solutions in the intervals 

(3.2) [x c + 6u c (t - t c ) - ae 6/7 , x c + 6u c (t - t c ) + ae 6/7 } 

where a is an e- independent constant (typically we take a = 3). 

Numerically we find at the critical time that the Loo-norm of the difference 
between Hopf and KdV solution scales like e a where a = 0.2869 (2/7 = 0.2857 . . .) 
with correlation coefficient r = 0.9995 and standard error <r a = 6.9* 10 -3 . Thus we 
confirm the expected scaling behavior within numerical accuracy. We also test this 
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difference for times close to t c . The relations (12.11|) suggest, however, a rescaling 
of the time, i.e., to compare solutions for different values of e at the same value 
of T. We compute the respective solutions for KdV times t±(e) = t c ± O.le 4 / 7 . 
Before breakup at t_ we obtain a — 0.31 with r — 0.999 and a a = 9.8 * 10~ 3 , i.e., as 
expected a slightly larger value than 2/7. After breakup at t+ we find a = 0.26 with 
0.9995 and a a — 6.6* 10~ 3 . We remark that after the breakup time, the asymptotic 
solution is obtained by gluing the Hopf solution and the theta-functional solution 

(G3). 

These results indicate that the scalings in (|2.1ip are indeed observed by the 
KdV solution. We show the corresponding situation for t- for two values of e in 
Fig. [3] 



e=10" 1 

0i 1 1 1 1 r- 




-| I i i i i i i i i 

-1.524 -1.522 -1.52 -1.518 -1.516 -1.514 -1.512 -1.51 

x 



FIGURE 3. KdV solution (blue) and Hopf solution (green) at the 
times t—(e) in a rescaled interval for two values of e. 

3.2. Multiscale solution. In Fig. [4] we show the numerical solution of the 
KdV equation for the initial data uq and the corresponding PI2 solution (|2.11[) 
for e = 10~ 2 close to breakup. It can be seen that the PI2 solution (|2.1ip gives 
a correct description of the KdV solution close to the breakup point. For larger 
values of \x — x c \ the multiscale solution is not a good approximation of the KdV 
solution. A similar situation is shown in Fig. [5] for the case e = 10~ 3 . Obviously 
the approximation is better for smaller e. Notice that the asymptotic description 
is always better near the leading edge than near the trailing edge. 

In Fig. [6] we plot in green the difference between the PI2 multiscale solution 
and the KdV solution and in blue the difference between the KdV solution and the 
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Figure 4. The blue line is the solution of the KdV equation for the 
initial data uq(x) — —1/ cosh 2 x and e = 10 -2 , and the green line 
is the corresponding multiscale solution given by formula (|2.11|) . 
The plots are given for different times near the point of gradient 
catastrophe (x c ,t c ) of the Hopf solution. Here x c ~ —1.524, t c ~ 
0.216. 

asymptotic solutions (|2.3p and (|2.4p . It is thus possible to identify a zone around x c 
in which the multiscale solution gives a better asymptotic description. The limiting 
values of this zone rescaled by x c are shown in Fig. [7] for the critical time. It can 
be seen that the zone always extends much further to the left (the direction of 
propagation) than to the right. The width of this zone scales roughly as e 3 / 7 , more 
precisely we find e a with a — 0.468, r = 0.981 and <r a = 0.073. We observe that 
the numerical scaling is smaller than the one predicted by the formula (|2 . 1 2|) . The 
matching of the multiscale and the Hopf solution can be seen in Fig. [H 

For larger times, the asymptotic solution (|2.3|) and (|2.4p gives as expected a 
better description of the KdV solution, see Fig. [5] for e = 10~ 3 and t = 0.226. 
Close to the leading edge, the oscillations are, however, better approximated by the 
multiscale solution. 

To study the scaling of the difference between the KdV and the multiscale 
solution, wc compute the norm of the difference between the solutions in the 
rescaled x-interval (|3.2|) with a = 3. Wc find that this error scales at the critical 
time roughly like e 5 / 7 . More precisely we find a scaling e a where a = 0.708 (5/7 = 
0.7143 . . .) with correlation coefficient r = 0.9998 and standard error a a — 0.012. 
Before breakup at the times t-(e) we obtain a = 0.748 with r = 0.9996 and a a = 
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Figure 5 . The blue line is the solution of the KdV equation for the 
initial data uq(x) — —1/ cosh 2 x and e = 10 -3 , and the green line 
is the corresponding multiscale solution given by formula (|2.11|) . 
The plots are given for different times near the point of gradient 
catastrophe (x c ,t c ) of the Hopf solution. 



0.016, after breakup at the times t+(e) we get a — 0.712 with r = 0.9999 and 
er a = 6.2 * 10~ 3 . Notice that the values for the scaling parameters are roughly 
independent of the precise value of the constant a which defines the length of the 
interval (|3.2p . For instance for a = 2, we find within the observed accuracy the same 
value. In [4] Claeys and Vanlessen showed that the corrections to the multiscale 
solution appear in order e 3 / 7 . For the values of e we could study for our KdV 
example, the corrections are apparently of order e 5 / 7 . 



4. Outlook 

The Camassa-Holm equation [3] (see also [12j ) 

(4.1) u t + 6uu x - e 2 {u xxt + 4u x ii xx + 2uu xxx ) = 

admits a bi-Hamiltonian description after the following Miura-type transformation 

(4.2) m = u — e 2 u xx . 
One of the Hamiltonian structure takes the form 

(4.3) {m(x), m(y)} = S'(x - y) - e 2 S"'(x - y) 
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Figure 6. The blue line is the difference between the solution of 
the KdV equation for the initial data uq(x) = —1/ cosh 2 x and 
e = 10~ 2 and the multiscale solution, and the green line is the 
difference between the asymptotic solution and the KdV solution. 
The plots are given for different times near the point of gradient 
catastrophe (x c , t c ) of the Hopf solution. 



so that the Camassa-Holm flow can be written in the form 
(4.4) m t = {m(x),H}, H = j ' (u 3 + uu 2 x )dx. 

To compare the Hamiltonian flow in (|1.2p with the one given in (|4.4p one must first 
reduce the Poisson bracket to the standard form {u(x),u(y)}i = S'(x — y) by the 
transformation 

u = (1 - e d x ) m = m+ -e m xx + -e m xxxx + .. .. 

Z a 

After this transformation, the Camassa-Holm equation will take for terms up to 
order e 4 the form 

~h 6ti Ux ~\~ £ ^U X U XX ~\- All U xxx 'j € (20 Uxx^xxx "f" 12 Ux^xxxx ~l _ 4li Uxxxxx) ' ' 0. 

which is equivalent to (|1.2p after the substitution 
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Figure 7. Limiting values of the zone where the multiscale solu- 
tion provides a better asymptotic description of the KdV solution 
than the Hopf solution for t = t c . The x values are rescaled with 




Figure 8. Difference of the KdV and the multiscale solution 
(blue) and the KdV and the Hopf solution (green) for the initial 
data uq(x) = —1/ cosh x at t = t c for two values of e. 



At the critical point (x c ,t c , u c ) the Camassa-Holm solution behaves according to 
the conjecture in |10j as 

u(x,t,e) = u c - (^f) ' U{X,T) + 0{^), c = 4u c 
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Figure 9. The blue line is the difference between solution of the 
KdV equation for the initial data uq(x) — —1/ cosh 2 a; and e = 
10 -3 and the multiscale solution, and the green line is the difference 
between the asymptotic solution and the KdV solution. The plots 
are given for t = 0.226. 

where 

1 / e \ 1/7 ( 1 \ 1/7 

x= -e{m) ^-^-^ -*-))« T= {7wj {t - tc) 

In Fig. [TD] we show the numerical solution to the CH equation for the initial 
data uo = — scch 2 (a;) and e = 1CP 2 at several values of time near the point of 
gradient catastrophe of the Hopf equation. It is interesting to compare this to the 
corresponding situation for the KdV equation in Fig. [4] It can be seen that there 
are no oscillations of the CH equation on left side (the direction of the propagation) 
of the critical point, whereas in the KdV case all oscillations are on this side. The 
quality of the approximation of the CH and the KdV solution by the multiscale 
solution is also different. In the KdV case, the solution is well described by the 
multiscale solution on the leading part which includes the oscillations, whereas the 
approximation is less satisfactory on the trailing side. A similar behavior is observed 
in the CH case, but since the oscillations are now on the trailing side, they are not 
as well approximated as in the KdV case. The leading part of the solution near the 
critical point is, however, described in a better way. 
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Figure 10. The blue line is the solution of the CH equation for 
the initial data Uq(x) = — 1/cosh 2 x and e = 10~ 2 , and the green 
line is the corresponding multiscale solution. The plots are given 
for different times near the point of gradient catastrophe (x c , t c ) of 
the Hopf solution. Here x c ~ —1.524, t c ~ 0.216. 

The same qualitative behavior can also be observed for smaller e in Fig. QTJ 
though the quality of the approximation increases as expected on the respective 
scales. Note that we plotted in Fig. [TUland Fig. [TT]the CH solution instead of the 
function u, since there are no visible differences between the two for the used values 
of e. 

Appendix A. Numerical solution of the fourth order equation 

We are interested in the numerical solution of the fourth order ordinary equa- 
tion (ODE) l|1.5p with the asymptotic conditions (|1.6[) . Numerically we will consider 
the equation on the finite interval [Xj,X r ], typically X r = —Xi = 100. In the ex- 
terior of this interval the solution to the equation Q1.5P is obtained in the form of a 
Laurent expansion of F around infinity in terms of Y = X 1 ^ , 

(A.1) U = Y + ±t^. 

n=l 

We find the non-zero coefficients (not-given coefficients vanish) a\ = 2T, a§ = 
-8T 3 /3, a 6 = 1/18, a 7 = 16T 4 /3, a 8 = -5T/27, a w = 14T 2 /27, a n = -256T 6 /9, 
012 = 16T 3 /3, a 13 = 640T 7 /9 - 7/108, ...This expansion also determines the 
boundary values we impose at Xi, X r for U and Ux- 

The solution in the interval [Xi , X r ] is numerically obtained with a finite differ- 
ence code based on a collocation method. The code bvp4c distributed with Matlab, 
see |19j for details, uses cubic polynomials in between the collocation points. The 
ODE (|1.5p is rewritten in the form of a first order system. With some initial guess 
(we use Uq = —X 1 / 3 as the initial guess), the differential equation is solved itera- 
tively by linearization. The collocation points (we use up to 10000) are dynamically 
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Figure 11. The blue line is the solution of the CH equation for 
the initial data Uq(x) = — 1/cosh 2 x and e = 10~ 3 , and the green 
line is the corresponding multiscale solution. The plots are given 
for different times near the point of gradient catastrophe (x c , t c ) of 
the Hopf solution. 



adjusted during the iteration. The iteration is stopped when the equation is sat- 
isfied at the collocation points with a prescribed relative accuracy, typically 10~ 6 . 
The values of U in between the collocation points are obtained via the cubic poly- 
nomials in terms of which the solution has been constructed. This interpolation 
leads to a loss in accuracy of roughly one order of magnitude with respect to the 
precision at the collocation points. To test this we determine the numerical solution 
via bvp4c for (|1.5|) on Chebychev collocation points and check the accuracy with 
which (|1.5p is satisfied via Chebychev differentiation, see e.g. [22] . We are inter- 
ested here in values of \T\ < 1 and \X\ < 10. It is found that the numerical solution 
with a relative tolerance of 10 -6 on the collocation points satisfies the ODE to the 
order of better than 10 -4 , see Fig. IT21 where we show the residual A by plugging the 
numerical solution into the differential equation. It is straight forward to obtain 
higher accuracy by requiring a lower value for the relative tolerance, but we will 
only need an accuracy of the solution of the order of 10 -4 here. 
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